Tutorial Part IV - Movement - Solution Exercise 2
Exercise 4.2
Three animals have been collared and tracked, and we want to know their home range size and whether the animals actually avoid agriculture and prefer forest.
There is telemetry data available from three animals (DieTiere_32633.shp) in the folder data-raw. Additionally, land cover (Landnutzung_3035.shp) and rivers (linearegew_32633.shp) are available in the geo-raw folder. Land use has a column ‘NS1’ with coded info on 1 (urban areas), 2 (agriculture), 3(forest), 4(swamps), 5 (water bodies).
After loading libraries, the data and making a plot, please do the following:
a) Calculate their home ranges.
b) Buffer the home range with the mean displacement per step of all three animals
c) Clip the home ranges with the underlying landscape and calculate the percentage
of land use classes 2/ agriculture and 3/ forest within the home range (hint: use st_intersection)
d) In which land use are the animals actually located? Extract land use underneath each telemetry location (hint: st_join). Do the animals actually avoid agriculture?
Get started
# Open libraries
## general basics
library(here)
library(ggplot2)
## the basic packages for spatial analyses and viz
library(sf)
library(terra)
library(tmap)
## for telemetry and home range analyses
library(move)
library(adehabitatHR)
library(sp)Define workspace and folder structure
work_wd <- here()
data_raw_wd <- here('data-raw')
geo_raw_wd <- here('data-raw/geo-raw')
output_wd <- here('output')Load data
# if you are not sure about the file names, check the folder content with
# list.files(data_raw_wd)
animals <- st_read(dsn=data_raw_wd,layer="DieTiere_32633") ## Reading layer `DieTiere_32633' from data source
## `C:\Users\kramer\PopDynIZW Dropbox\Steph Kramer\_GitHub\Course4_MoveQ\data-raw'
## using driver `ESRI Shapefile'
## Simple feature collection with 388 features and 2 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 390397.5 ymin: 5862042 xmax: 396748.2 ymax: 5868366
## Projected CRS: WGS 84 / UTM zone 33N
landcov <- st_read(layer="Landnutzung_3035", dsn=geo_raw_wd) ## Reading layer `Landnutzung_3035' from data source
## `C:\Users\kramer\PopDynIZW Dropbox\Steph Kramer\_GitHub\Course4_MoveQ\data-raw\geo-raw'
## using driver `ESRI Shapefile'
## Simple feature collection with 182 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 4539362 ymin: 3304861 xmax: 4566087 ymax: 3324860
## Projected CRS: ETRS89-extended / LAEA Europe
rivers <- st_read(layer="linearegew_32633", dsn=geo_raw_wd) ## Reading layer `linearegew_32633' from data source
## `C:\Users\kramer\PopDynIZW Dropbox\Steph Kramer\_GitHub\Course4_MoveQ\data-raw\geo-raw'
## using driver `ESRI Shapefile'
## Simple feature collection with 2834 features and 13 fields
## Geometry type: MULTILINESTRING
## Dimension: XYZ
## Bounding box: xmin: 381989.2 ymin: 5852103 xmax: 408856.6 ymax: 5871804
## z_range: zmin: 0 zmax: 0
## Projected CRS: WGS 84 / UTM zone 33N
Define or transform CRS
lc_32633 <- st_transform(landcov, 32633)Make a plot of the data
#Create a color code for the 5 land use classes
mycol <- c('red','antiquewhite','chartreuse4','cornsilk4','blue')
## TODO - argh, Punkte und Linien erscheinen nicht im html....
plot(lc_32633[,'NS1'], col=mycol[lc_32633$NS1],border='transparent')
plot(animals[,1],col=animals$ID,add=TRUE)
plot(rivers[,1], col= 'blue', lwd=0.1,add=T)#Create a color code for the 5 land use classes
mycol <- c('red','antiquewhite','chartreuse4','cornsilk4','blue')
# TODO - plot mit colorcode mycol funktioniert nicht....
ggplot(lc_32633) +
geom_sf(data=lc_32633[,'NS1'], aes(fill = mycol[NS1]), color='transparent') +
geom_sf(data=animals) +
geom_sf(data=rivers, color='blue') +
theme_classic()#Create a color code for the 5 land use classes
mycol <- c('red','antiquewhite','chartreuse4','cornsilk4','blue')
myid <- c('yellow','pink','red')
tmap_mode(mode = "view")
tm_shape(shp = lc_32633[,'NS1']) +
tm_polygons(col = "NS1", palette = mycol) +
tm_shape(shp = animals[,'ID']) +
tm_dots(col = "ID", palette = myid) +
tm_shape(shp = rivers[,1]) +
tm_lines(col = 'blue') a) & b)
First , we retrieve the steplength. For this, transform the sf-Object into a move-Object.
anim_coords <- as.data.frame(st_coordinates(animals)) #returns a list
animals$X <- anim_coords$X
animals$Y <- anim_coords$Y
animals_ng <- st_drop_geometry(animals)
animals_move <- move(x=animals_ng$X, y=animals_ng$Y,
time = as.POSIXct(animals_ng$MYDATUM,
format = "%Y-%m-%d 8:0:0", tz="CET"),
proj=CRS("+init=epsg:32633"),
data=animals_ng, animal = animals_ng$ID, sensor='GPS')
plot(animals_move, type='l') # since we have 3 animals, a list with 3 elements is returned
steplength_per_animal <- distance(animals_move)
mean_stepl_p_anim <- lapply(steplength_per_animal,mean)
(mean(unlist(mean_stepl_p_anim)))## [1] 303.0478
The mean step length is approx. 300 m. Now use this value to buffer the MCPs.
# for using mcp() in {adehabitatHR}, convert to SpatialPointsDataframe of {sp}
animals_sp <-as(animals, "Spatial")
mcp <- mcp(animals_sp[,'ID'], percent=95)
# the buffer command needs however an sf object
mcp_sf <- as(mcp, "sf")
mcp_buf <- st_buffer(mcp_sf, 300)
# the plot is correct
plot(mcp)
plot(animals_sp, add=TRUE, col= animals_sp$ID)
plot(mcp_buf[,1], add=TRUE)
plot(animals_sp, add=TRUE, col= animals_sp$ID)
plot(mcp, add=TRUE)# save buffer
st_write(obj = mcp_buf,
dsn = output_wd,
layer = 'mcp_buffer_32633',
driver = 'ESRI Shapefile',
delete_layer = TRUE)## Deleting layer `mcp_buffer_32633' using driver `ESRI Shapefile'
## Writing layer `mcp_buffer_32633' to data source
## `C:/Users/kramer/PopDynIZW Dropbox/Steph Kramer/_GitHub/Course4_MoveQ/output' using driver `ESRI Shapefile'
## Writing 3 features with 2 fields and geometry type Polygon.
c)
# check the total area of the buffer in km2 per animal:
buf_area <- st_area(mcp_buf)/1e6
# clip with land use
mcp_buf_lc <- st_intersection(lc_32633, mcp_buf)
# have a look at the object:
mcp_buf_lc## Simple feature collection with 45 features and 14 fields
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: 390158.2 ymin: 5861742 xmax: 397048.2 ymax: 5868666
## Projected CRS: WGS 84 / UTM zone 33N
## First 10 features:
## AREA PERIMETER CLCGES_BES CLCGES_B_1 NS BRD SZENE SPHEROID TKNR NS1
## 2 33397200 79888.10 23163 23162 312 1 19323 Krassovsky N33111 3
## 5 11562800 44001.70 23265 23264 211 1 19323 Krassovsky N33111 2
## 15 18208500 73976.60 23764 23763 231 1 19323 Krassovsky N33111 2
## 16 8568800 37160.80 23788 23787 231 1 19323 Krassovsky N33111 2
## 28 5142930 13701.10 24028 24027 312 1 19323 Krassovsky N33111 3
## 35 1227710 4956.00 24314 24313 243 1 19323 Krassovsky N33111 2
## 36 1448920 11624.70 24348 24347 311 1 19323 Krassovsky N33111 3
## 41 654073 3560.28 24431 24430 211 1 19323 Krassovsky N33111 2
## 42 3534710 11425.00 24454 24453 311 1 19323 Krassovsky N33111 3
## 44 46474200 125933.00 24480 24479 211 1 19323 Krassovsky N33111 2
## NS2 NS3 id area geometry
## 2 31 312 1 1444.323 MULTIPOLYGON (((394196.9 58...
## 5 21 211 1 1444.323 POLYGON ((392577.5 5866547,...
## 15 23 231 1 1444.323 POLYGON ((391413.6 5866716,...
## 16 23 231 1 1444.323 POLYGON ((395139.3 5866196,...
## 28 31 312 1 1444.323 POLYGON ((390833 5866283, 3...
## 35 24 243 1 1444.323 POLYGON ((392292.4 5866469,...
## 36 31 311 1 1444.323 POLYGON ((395809.1 5865790,...
## 41 21 211 1 1444.323 POLYGON ((394729.2 5865833,...
## 42 31 311 1 1444.323 POLYGON ((390395.2 5865319,...
## 44 21 211 1 1444.323 MULTIPOLYGON (((392247.2 58...
# take care, you have to recalculate the area!
mcp_buf_lc$areakm2 <- st_area(mcp_buf_lc)/1e6
# now, for each animal id, summarize the area per NS1-category 2 and 3:
agri <- tapply(mcp_buf_lc$areakm2[mcp_buf_lc$NS1 == 2],
mcp_buf_lc$id[mcp_buf_lc$NS1 == 2],sum)
forest <- tapply(mcp_buf_lc$areakm2[mcp_buf_lc$NS1 == 3],
mcp_buf_lc$id[mcp_buf_lc$NS1 == 3],sum)
# percentages
round((agri/buf_area*100),0)## Units: [1/m^2]
## 1 2 3
## 71 48 85
round((forest/buf_area*100),0)## Units: [1/m^2]
## 1 2 3
## 29 52 14
d)
spatjoin <- st_join(animals,lc_32633)
myres <- table(spatjoin$NS1,spatjoin$ID)
myres##
## 1 2 3
## 2 48 25 95
## 3 147 64 9
# now divide that by the total number of locations per animal
myn <- table(animals$ID)
round((myres[1,]/myn*100),0) # agriculture##
## 1 2 3
## 25 28 91
round((myres[2,]/myn*100),0) # forest##
## 1 2 3
## 75 72 9
# compare the values with the percentages aboveFor example, animal 1 had 71% of agriculture and 29% of forest in its HR, however, only 25% of the locations were actually in agriculture, but 75% of its locations were in forest. We could conclude that animal 1 avoids agriculture.
END
Session Info
Sys.time()## [1] "2022-03-07 02:33:24 CET"
sessionInfo()## R version 4.1.2 (2021-11-01)
## Platform: i386-w64-mingw32/i386 (32-bit)
## Running under: Windows 10 x64 (build 19043)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_Germany.1252 LC_CTYPE=English_Germany.1252
## [3] LC_MONETARY=English_Germany.1252 LC_NUMERIC=C
## [5] LC_TIME=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] adehabitatHR_0.4.19 adehabitatLT_0.3.25 CircStats_0.2-6
## [4] boot_1.3-28 MASS_7.3-54 adehabitatMA_0.3.14
## [7] ade4_1.7-18 deldir_1.0-6 move_4.1.6
## [10] rgdal_1.5-28 raster_3.5-11 sp_1.4-6
## [13] geosphere_1.5-14 tmap_3.3-2 terra_1.4-22
## [16] sf_1.0-5 ggplot2_3.3.5 here_1.0.1
##
## loaded via a namespace (and not attached):
## [1] RColorBrewer_1.1-2 httr_1.4.2 rprojroot_2.0.2
## [4] tools_4.1.2 bslib_0.3.1 utf8_1.2.2
## [7] R6_2.5.1 KernSmooth_2.23-20 DBI_1.1.1
## [10] colorspace_2.0-2 withr_2.4.3 tidyselect_1.1.1
## [13] leaflet_2.0.4.1 compiler_4.1.2 leafem_0.1.6
## [16] textshaping_0.3.6 xml2_1.3.3 bookdown_0.24
## [19] sass_0.4.0 scales_1.1.1 classInt_0.4-3
## [22] proxy_0.4-26 systemfonts_1.0.3 stringr_1.4.0
## [25] digest_0.6.29 rmarkdown_2.11 base64enc_0.1-3
## [28] dichromat_2.0-0 pkgconfig_2.0.3 htmltools_0.5.2
## [31] highr_0.9 fastmap_1.1.0 htmlwidgets_1.5.4
## [34] rlang_0.4.12 farver_2.1.0 jquerylib_0.1.4
## [37] generics_0.1.1 jsonlite_1.7.3 crosstalk_1.2.0
## [40] dplyr_1.0.7 magrittr_2.0.1 s2_1.0.7
## [43] Rcpp_1.0.8 munsell_0.5.0 fansi_1.0.2
## [46] abind_1.4-5 lifecycle_1.0.1 stringi_1.7.6
## [49] leafsync_0.1.0 yaml_2.2.1 tmaptools_3.1-1
## [52] grid_4.1.2 parallel_4.1.2 crayon_1.4.2
## [55] lattice_0.20-45 stars_0.5-5 knitr_1.36
## [58] pillar_1.6.4 codetools_0.2-18 wk_0.5.0
## [61] XML_3.99-0.8 glue_1.6.0 evaluate_0.14
## [64] leaflet.providers_1.9.0 png_0.1-7 vctrs_0.3.8
## [67] rmdformats_1.0.3 gtable_0.3.0 purrr_0.3.4
## [70] assertthat_0.2.1 cachem_1.0.6 xfun_0.28
## [73] lwgeom_0.2-8 e1071_1.7-9 ragg_1.2.1
## [76] class_7.3-19 viridisLite_0.4.0 tibble_3.1.6
## [79] memoise_2.0.1 units_0.7-2 ellipsis_0.3.2